library(desolve)

Stability<-function(t,state,parameters){
  C1=state[1]
  C2=state[2]
  C3=state[3]
  k1=exp(parameters[1])
  k2=exp(parameters[2])
  dC1<--k1*C1*C3
  dC2<- k2*C3
  dC3<- k1*C1*C3-k2*C3
  list(c(dC1, dC2, dC3))
}

likelihood<-function(x, ini, dataX, dataY){
  n=ncol(dataX)
  prob_C2=0
  for(i in 1:n){
    y=c(ini[i,],100-sum(ini[i,]))
    times=c(0,dataX[,i])
    pred=ode(y,times, func=Stability,parms=x) 
    C2_pred=pred[-1,3]
    residual=log(dataY[,i])-log(C2_pred)
    prob_C2=prob_C2+sum(dnorm(as.matrix(residual), mean=0, sd=exp(x[length(x)]), 
      log=TRUE))
   }      
   prob_C2
}


system.time(
  post.samp <- MCMCmetrop1R(likelihood, 
    theta.init=c(-9.29,-4.09,-1.8959),
    thin=1, mcmc=1e4, burnin=0.25e4, tune=rep(1.5,3), verbose=10,
    logfun = TRUE, force.samp = TRUE, V = NULL, optim.method = "SANN", 
    ini=initial_conditions, dataX=dataX, dataY=dataY)  
, gcFirst = TRUE) 


likelihood <- function (x, ini, dataX, dataY1, U = NULL, center = NULL, 
	pred = FALSE, uniforming = 0, raw_val = TRUE)
{
	if (is.null (U)) {
    if(is.null(center)){
		  xt <- rep (0, length (x))
		  #PARAMETERS
		  xt <- x}else{
      xt <- x+center}
	} else {
    if(is.null(center)){                         
      xt <- Solve (U, x)
    }else{
		  xt <- Solve (U, x)
    xt <- xt + center  }
	}
  n=ncol(dataX)
  prob_C2=0
  for(i in 1:n){
    y=c(ini[i,],100-sum(ini[i,]))
    times=c(0,dataX[,i])
    pred_out=ode(y, times, func=Stability, parms=xt) 
    C2_pred=pred_out[-1,3]
    residual=log(dataY1[,i])-log(C2_pred)
    prob_C2=prob_C2+sum(dnorm(as.matrix(residual), mean=0, sd=exp(x[length(x)]), 
      log=TRUE))
  }
	if (raw_val == TRUE) {
		if (pred == FALSE) {
			c (xt, (uniforming + prob_C2))
		} else {
      		list(prob_C2, y1_pred, xt, exp(x[length(x)]))
		}               
	} else {
		if (pred == FALSE) {
			c (x, (uniforming + prob_C2))
		} else {
      		list (prob_C2, y1_pred, x, exp(x[length(x)]))
		}
	}
}
    

    
  